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Abstract 



We study the problem of optimal pricing and hedging of a European option written 
on an illiquid asset Z using a set of proxies: a liquid asset S, and N liquid European 
options Pi, each written on a liquid asset Yi,i = 1,N. We assume that the S'-hedge 
is dynamic while the multi-name Y-hedge is static. Using the indifference pricing 
I/"") ■ approach with an exponential utility, we derive a HJB equation for the value function, 

and build an efficient numerical algorithm. The latter is based on several changes 
of variables, a splitting scheme, and a set of Fast Gauss Transforms (FGT), which 
turns out to be more efficient in terms of complexity and lower local space error than a 
finite-difference method. While in this paper we apply our framework to an incomplete 
market version of the credit-equity Merton's model, the same approach can be used for 
other asset classes (equity, commodity, FX, etc.), e.g. for pricing and hedging options 
with illiquid strikes or illiquid exotic options. 



1 Introduction 

This work is an extension of our paper j3| where the following problem was considered. 
We support a trader who wants to buy (or sell) a European option Cz on asset Z with 
maturity T and payoff Gz- The trader wants to hedge this position, but the underlying 
asset Z is illiquid. However, some liquid proxies of Z are available in the marketplace. First, 



* Opinions expressed in this paper are those of the authors, and do not necessarily reflect the view of 
JPMorgan Chase and Numerix 
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there is a financial index (or simply an index) S (such as e.g. S&P500 or CDX.NA)Q| whose 
market price is correlated with Z. In addition, there is another correlated asset Y which 
has a liquidly traded option Cy with a payoff Cy similar to that of Cz, and with the same 
maturity T. The market price py of Cy is also known. 

Our trader realizes that hedging Z-derivative with the index S alone may not be suffi- 
cient for a number of reasons. First, she might be faced with a situation where correlation 
coefficients p yz , p sz (which for simplicity are assumed to be constant) are such that p yz > p sz . 
In this case we would intuitively expect a better hedge produced by using Y or Cy as the 
hedging instruments. Second, if we bear in mind a stochastic volatility-type dynamics for 
Z, the stochastic volatility process may be "unspanned", i.e. the volatility risk of the option 
may not be traded away by hedging in option's underlying^]. If that is the case, one might 
want to hedge the unspanned stochastic volatility by trading in a "similar" option on the 
proxy asset Y. So our trader is contemplating a hedging strategy that would use both S and 
Y. To capture an "unspanned" stochastic volatility, the trader wants to use a derivative Cy 
written on Y rather than asset Y directly. 

As transaction costs are usually substantially higher for options than for underlyings, 
our trader sets up a static hedge in Cy and a dynamic hedge in S t . The static hedging 
strategy amounts to selling (or buying) a units of Cy options at time t = 0. An optimal 
hedging strategy would be composed of a pair (a*, 7r*) where a* is the optimal static hedge, 
and 7r* (where < s < T) is an optimal dynamic hedging strategy in index St. The pair 
(a*, Ti*) should be obtained using a proper model. The same model should produce the 
highest/lowest price for which the trader should agree to buy or sell the Z-option. 



In [7| we developed a model that formalizes the above scenario by supplementing it with 
the specific dynamics for asset prices St, Y t and Z t , and providing criteria of optimality for 
pricing options Cz- For the former, we use a standard correlated log- normal dynamics. For 
the latter, we employ the utility indifference framework with an exponential utility, pioneered 
by @, EH and others, see e.g. j^f for a review. We showed that this results into a tractable 
setup with analytical (in quadratures) expressions for optimal hedges and option prices. For 
more details and links to the related literature, see Q. 

To extend this model, we notice that availability of just one asset for the static hedge in 
our model is very restrictive. More generally, we may assume that N liquid options on assets 
Yi, i = 1, . . . , N are available in the marketplace, where all Y$ have similar correlations p ZiVi 
with asset Z t . Therefore, all N options could be used in this scenario to set up for static 
hedging of the Z-option. All in total, we have N + 1 assets for a static-dynamic hedge 
optimization problem. 

This is the problem addressed by the present work. Similar to [7] , we use the indifference 
pricing approach and an exponential utility function to derive a HJB equation for the value 
function. In the present case, the HJB equation is (N + l)-dimensional. We develop an 
efficient numerical algorithm to solve the HJB equation. Our approach is based on several 

1 Here we refer to this instrument as an index, but it could be any "linear" instrument such as stock, 
forward, etc. 

2 The notion of unspanned stochastic volatility was introduced in jij. For a discussion of such scenarios 
for e.g. commodities markets, see [5oj |. 
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changes of variables, a splitting scheme, and a set of Fast Gauss Transforms (FGT), which 
turns out to be more efficient in terms of complexity and lower local space error than a 
finite-difference method. 

Before presenting our notation and convention, we note that while the mathematic frame- 
work developed below is general and can be applied to various asset classes, for definiteness 
below we follow Ref. [7J and specialize on pricing and hedging on illiquid debt within a 
version of the Merton credit-equity model. This setting might be of interest for modeling 
counterparty value adjustments (CVA) for over-the-counter (OTC) derivatives. However, a 
similar framework can be developed for other cases where the utility-indifference approach 
is useful, e.g. executive stock options could be priced along the same lines. 

We therefore assume a market where the following instruments can be traded: 

• a risk-free zero-coupon bond Bq 

• a risky non- default able index S 

• a set of liquid bonds By. issued by firms Yi, % = 1, N, with market prices p Y . 

Our financial model for pricing and hedging an illiquid credit Bz amounts to computing its 
price pz in terms of all py t and S at t — 0, along with optimal hedges. Note that as long as 
issuers of Y$ and Z are imperfectly correlated, the liquid bonds By. provide only a partial 
hedge for Bz- 

As we are in the incomplete market setting, risk of Z cannot be perfectly hedged by 
(By t ,S), hence both the price and hedge ratios will be different for different investors, de- 
pending on their risk preferences and a (non-unique) hedging strategy. 

Therefore, the idea is to hedge an exposure to a counterparty with illiquid credit (a long 
position in bond B z ) by taking static short positions in a set of proxy liquid debts By., plus 
possibly using a dynamic trading strategy 6 t in the index S. 

Assume we statically hedge bond Bz by selling a, zero-coupon bonds issued by firm Yi 
for their market price py v The cash amount available for investing in bonds and index is 
x + dipy^ where x is the initial cash minus the price paid for B Z - 

We further use an indifference pricing principle to derive an HJB equation which describes 
an evolution of the investor utility function in our setup. As no closed form solutions are 
known for the utility indifference pricing with N > 2, we suggest a very efficient numerical 
method of solving the HJB using a combination of a special change of variables and a 
particular splitting scheme. Its total complexity is 0(M + N)M + 0(M + 1) ps M 2 , M > N, 
where M is the number of nodes for Fast Gauss Transform (FGT) used in calculations. This 
is significantly less than e.g. the total complexity 0(N M ) of finite difference methods. Both 
the theoretical setup and numerical algorithms presented below are the main results of this 
paper, which to our knowledge are new. 

The rest of the paper is organized as follows. Section [2] describes a general setup of the 
problem and indifference pricing framework. In section [3] we derive a corresponding HJB 
equation for our model. Section H] introduces new factorized (adiabatic) variables, and shows 
that in new variables the HJB equation transforms to a iV-dimensional heat equation with 
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an extra non-linear term. This term is proportional to $^ o /$, i.e. it contains only the first 
derivative of the dependent variable $ wrt the first independent variable yo but no other 
derivatives. We describe an efficient numerical algorithm to compute coefficients of such a 
transformation. The next section shows how the transformed HJB equation can be solved 
numerically using Strang's splitting. We show that the problem reduces to the solution of one 
iV-dimensional and two one-dimensional heat equations. For the latter task, we show how 
to use FGT to decrease the total complexity of the method. Section [H] discusses calibration 
of the method to the market data. The final section concludes. 



2 Static hedging in indifference pricing framework 



Borrowing from an approach of [12[ for a similar (but not identical) setting, we now show 
how the method of indifference utility pricing can be generalized to incorporate our scenario 
of a mixed dynamic-static hedge. 

To this end, let U.(Yr, Zt) be the final payoff of the portfolio consisting of our option 
positions, i.e. 

n a (Y T ,Z T ) = G z -J2®iG Yi (1) 

i 

For convenience let us further denote asset Z as Yq. As long as all European options 
Cy v i G [0, N] pay at the same maturity T, we can view this as the payoff of a combined 
("static hedge portfolio") option g(a , a^), which involves payoffs Gy i of all derivatives 
Cy i . Such option may be priced using the standard utility indifference principle. The latter 
states that the derivative price g(ao, ...,a7v) is such that the investor should be indifferent 
to the choice between two investment strategies. With the first strategy, the investor adds 
the derivatives to her portfolio of bonds and stocks (or indices^) S, thus taking g(a , ajy) 
from, and adding ^\ «ipy 7 to her initial cash x. With the second strategy, the investor stays 
with the optimal portfolio containing bonds and the stocks/indices. 

The value of each investment is measured in terms of the value function defined as the 
conditional expectation of utility U(Wt) of the terminal wealth Wt optimized over trading 
strategies. In this work, we use an exponential utility function 



U(W) = -e 



-■yW 



where 7 is a risk-aversion parameter. In our case, the terminal wealth is given by the 
following expression: 

W T = X T + Tl ao >-> aN (Y 0:T ,...,Y N:T ) 

with Xt be the total wealth at time T in bonds and index S. In turn, the value function 
reads 

V(t,x,y ,...y N ) = (3) 



sup E 

TTt&M 



U(X T + U ao -^(Y 0tT ,...,Y N , T )) 



X t = x, Y 0tt = y , Y N>t = yN 



The stock is equivalent to our index S in the setting of the Merton's optimal investment problem. 
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where M. is a set of admissible trading strategies that require holding of initial cash x. The 
expectation in the Eq.([n]) is taken under the "real- world" measure P. 

For a portfolio made exclusively of stocks/indices and bonds, the value function for the 
exponential utility is known from the classical Merton's work: 

V°(x,t) = -e^ xerT -^ T (4) 

where r = T — t, r is the risk free interest rate assumed to be constant, and rj s = (/x s — r)/a s 
is the stock Sharpe ratio. 

In our setting, in addition to bonds and stocks/indices, we want to long Cy option and 
short aij units of every Cy x option to statically hedge our Cy position, or, equivalently, buy 
the g(ao, a at) option. 

From the Eq.([3]) the value function in our problem of optimal investment in bonds, index 
and the composite option g(a , a^v) has the following form: 



V(t,x,y , ...y N ) = (5) 

X t = x, Y , t = 2/o, — , Y Ntt = y N 



sup E 

TTt&M 



_ e -y(x T +n a o--<*N(Y 0>T ,...,YN, T )) 



where Xt is a cash equivalent of the total wealth in bonds and the index at time T. We 
represent it in a form similar to Eq. (j5j): 



V(t,x,y ,...y N ) = -e-T* e "^(t, y , ...y N ) (6) 

where function $ will be calculated in the next sections. The indifference pricing equation 
reads 

V(t, x, y , ...y N ) = V° ^t, x + g(a , a N ) - ^ oap Y A 
Plugging this in Eq.flU) and Eq.([n]) and re-arranging terms, we obtain 

1 N 
g(a , a N ) = e rT log $(r, y , ...y N ) + ^ a iPY,i 

The highest price of the Yo-derivative is given by choosing the optimal static hedge given by 
the numbers ai, ajv of the Yj-derivatives, i.e. 



g{a , ...,a Nj 



I - 

-e rT log $K--^)( T) y Qi ...y N ) + a iPY>i (7) 



a , ...,a 



ff) = arg max g(a* , ...,a* N ) 



where we temporarily introduced a superscripts a, in $( a S'-' a iv) to emphasize that the value 
function depends on all through a terminal condition. 
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3 The HJB equation 



To use the Eq.(J7]) and thus be able to compute both the option price and optimal static 
hedge, we need to find the "reduced" value function $. To accomplish this goal below we 
first derive the Hamilton- Jacobi-Bellman (HJB) equation for our model, and in the next 
section show how to efficiently solve it numerically. 

Let 9 be the investment strategy in the index. Optimal dynamic strategy can be obtained 
by using a general HJB principle 

V t + sup £"V = (8) 

where C n is the Markov generator, and tt = TT t (x) is the dynamic strategy at time t which 
depends on the initial cash amount x. 

Further assume that all state variables Sf,Yi,i G [0, N] follow a geometric Brownian 
motion process with time-dependent drifts and volatilities <7j, i G (x,0, ...,N) 

dS t = p x (t)S t dt + a x S t dWi x) 

dYi = faYidt + OiYtdW^ , i G [0, N\ 



Also following [17| assume that a riskless bond B t = 1 with maturity T is available for 
trading, yielding a constant interest rate r. Since our trading strategy implies a static 
position in all derivatives and dynamic positions in the index, real trading occurs in the time 
horizon [t, T],0 <t <T, and only between the two traded assets, i.e., the riskless bond B t 
and the risky asset St- If our total wealth at time t is X t = x and we invest amount 7r of this 
wealth into the index and the rest in a risk-free bond, the stochastic differential equation for 
X t is obtained as follows: 

dX t = r(X t - tt) dt + —dS t = (rX t + 7ca x r] s ) dt + na x dWf> , r] x = — 

ot cr x 

Then C 77 reads 

N 



1 N 

C* = [rx + n(fi x - r)] V x + -c^7rV xx + S P j p x ^ l o x Oy{Ky i V x , 



2 

i=0 

N . N N 

where V(t, x, yo, ...j/at) is defined on the domain K.(i, x, yo, y^) '■ [0, T] x [0, oo) x [0, oo) x 
... x [0, oo). 

Since C n is a regular function of tt, sup n is achieved at 

~r*(rr,\ — VxVx + ^2j = Q PxyjCiyiVxyj 

71 y x > ~ y 

v XX 
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Plugging this into Eq.flS]), we obtain 



N N N 

V t + rxV x + ViViVy,. + P^^oViV^v, 



(9) 



i=Q i=0 j=0 

^iV _ _ . t 7" \2 



1 (VxVx ~l~ X^i=Q Pxyi^yiViV xyn 

2 VCr 







This is a nonlinear PDE with respect to the dependent variable V(t, x, yo, i/n) with stan- 
dard boundary conditions (see [17]), and the terminal condition determined by a choice of 
the writer's maximal expected utility (value function) of the terminal wealth Wt- 

Note that so far the derivation is valid for a generic utility function. To make further 
progress, we specialize to the case of exponential utility in Eq.fl2]). The latter choice gives 
rise to a natural dimension reduction of the HJB equation. Indeed, the ansatz 



V(t, x, ?/o, y N ) = - exp (-<yxe r{T r) ) G(r, zq, z N ) 



(10) 



with Zi = \og(yi/Ki), i e [0, N] is both consistent with terminal condition Eq.flS]) and, upon 
substitution in flU, leads to a PDE for function G which does not contain variable x: 



^ n . n N I f N \ 2 

G T = --^IxG + ftG Vi + - 22 ^2 P'j a ' :i ( ' n-n. ~ 2Q\Z2 P*vi a iGyi I ) 

i=0 i=0 j=0 \i=0 J 



:ir 



where ft = ft - - r} x p xyi Oi. 

Equation Eq. ffTTl) is defined on the domain M.(t, zq, zn) : [0, T] x [—00,00) x ... x 
[—00, 00). The initial condition for this equation is obtained from Eq.flS]). 

In what follows, we choose a specific payoff of the form Eq.flT]) with IT = min(Yj, iQ), i e 
[0, N] where Ki are strikes. Then the terminal condition for G(r, z , z N ) reads 



G(0,z o , ...,z N ) = exp 
where z~ = min(zj,0). 



N 



"7 



K e z ° -^aiKie z 



i=i 



(12) 



4 The HJB equation and factorized variables 

The Eq. fllip is a (N + l)-dimensional parabolic equation with a non-linear (quadratic) term. 
No closed form solution is available for this case. Note that when N — 1, the HJB equation 
can be solved using an asymptotic expansion proposed in pj. Another relevant reference is 
To| that studies a related problem of counterparty risk of derivatives in incomplete markets 
with one traded and multiple non-traded assets^. 

4 We note that our splitting method (see below) is different from that used by the authors of [loj . In 
addition, their method is of the first order in time, while our method is of second order in time 
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Furthermore, straightforward applications of common numerical methods such as e.g. fi- 
nite differences would likely be inefficient in our setting. Indeed, assume that we approximate 
the non-linear term explicitly, as this does not affect stability of the FD scheme. Eq. ffTT]) 
then transforms to a (iV + l)-dimensional linear parabolic equation with a source term, which 
would be computationally costly to solve. 

An alternative to this solution, yet straightforward numerical approach could be con- 



structed as follows. We first use splitting (see e.g. |l6|) that reduces the original (N + 1)- 
dimensional problem to a set of iV + 1 one- dimensional problems. Thus, if every one- 
dimensional grid contains M nodes, and since every one- dimensional problem has a tridi- 
agonal matrix, the total complexity of the method is 0(M(N + 1)). Next we use the Fast 



Gauss Transform [21] to solve the resulting one- dimensional problems. 

Below we show that this straightforward approach can be significantly improved by 
rewriting the Eq. ffTTj) in new "factorized" variables. The reason that we call these vari- 
able "factorized" will be clear below. 

First, make a change of the dependent variable G — > <£> as follows: 

G(r, z , z N ) = e-^ r $(r, z , z N ), (13) 

so the first term in the rhs of the Eq. lfTTT) drops off the equation for $. 

Our further idea is to build a map y = (uo.-.i/n) — > u = (uo...un) such that in new 
variables u, both the Hessian matrix and the quadratic term in the Eq. fllip become diagonal. 

To be more specific, let us first introduce some matrix notation. Let A be the Hessian 
matrix, i.e. A = \\pijCTicrj\\, i,j £ [0, N]. Let a be a vector a = {pxy^i), i G [0, N]. Let R be 
a transformation matrix, i.e. u = R T y. Then we want to find such R that obeys 

R T AR = A, aR = B, 

where A is some diagonal matrix, and B = (1,0...0). 

In other words to determine d,2,--,dN we have to solve a system of non-linear algebraic 
equations aR = B with n = 2, . . . , iV wrt d.2, d^, where matrix R is defined implicitly via 
the solution of the eigenvalues problem R~ l DAR = A. This can be easily implemented, e.g. 
in Matlab just in few lines of code. The algorithm is pretty fast and converges to e = 10 -15 
within 25 msec for iV = 4 at Intel i7-2720 QM CPU 2.20 Ghz. 

Since eigenvectors are defined up to scaling, we fix it by choosing value —1 in the right- 
bottom corner of matrix D instead of adding an extra unknown d\. Accordingly, we solve a 
system of N — 1 equations C{ = 0, i & [2, N], rather then N equations. This results in the 
fact, that the first element of vector B could be whatever it becomes, rather than just 1. 

Based on Proposition 14.11 we conclude that the above algorithm transforms the Hessian 
matrix to the diagonal form. At the same time the last step of the algorithm guarantees 
that in new variables the quadratic form in the nominator of the non-linear term in Eq. (TTT]) 
contains just one (namely, the first) term. That is exactly what we wanted to achieve by 
doing the proposed change of variables. 

Some comments on the above algorithm should be made. First, logically the more our 
proxy assets correlate with the illiquid asset the better we can price the illiquid asset deriva- 
tives. This means that matrix \p\ has all elements, say in a range 0.5 < \pij\ < 1. Under these 
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Algorithm 1 Algorithm of building matrix of transformations R. 



1. Take a diagonal matrix D 
determined. 





■ 


■ 








d 3 ■ 
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■ 


■ d N 
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■ 





where d 2 , ...dN are the unknowns to be 



2. Assign some initial values to d 2 , ...d^ and solve an eigenvalues problem R -1 DAR = A, 
where A is a diagonal matrix with eigenvalues at the diagonal. Then use the following 
proposition 

Proposition 4.1 If R is a matrix of eigenvectors of DA, e.g. R~ l DAR = A, then 
R T AR is a diagonal matrix. 

Proof 1 See AppendixlAl 

3. Compute a vector C = \aR — B\. If all C 2 , are less then the method tolerance e 
- we are done. Otherwise take the next guess on d 2 , ...d N and proceed until converge. 



matrix DA becomes stiff with a high conditional number. Therefore, an accurate computa- 
tion of eigenvalues and eigenfunctions of such a matrix requires high precision arithmetics. 
That means that at a 32-bit architecture the proposed algorithm could fail to converge to 
the true solution with the required accuracy (despite it converges to some solutions with a 
bigger error). Moving the algorithm to a 64-bit architecture significantly improved conver- 
gency but still could fail when | close to 1. Therefore, in this case special algorithms 

of computing eigenvectors for stiff matrices have to be applied. 

After the transformation matrix is found we finally use a change of independent variables 

u = R T y + rM, M = ^- J fa (k)dk, p, N (k)dkj , 



to obtain 



This can also be written as 



2 / , r*~wm 2 $ 

i=0 



N 



$ r = ^£,;$ (15) 



i=0 



1 d 2 1 $ 2 Id 2 

2 ^ 2 $ 2 dyf 
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Here pi, % G [0, N] are the diagonal elements of the diagonal matrix R T AR (which is the 
Hessian matrix in new coordinates u), and bo is the first element of vector B. 

It is seen that in new variables operators Ci, i G [1, N] are linear. In addition, all operators 
G [0, N] are independent. That is why we call these new variables u as factorized. 

Example. Consider N = 3 and the following parameters of the model: 



Pvt. 



1.0 


0.9 


0.6 


0.5 


0.9 


1.0 


0.75 


0.7 


0.6 


0.75 


1.0 


0.6 


0.5 


0.7 


0.6 


1.0 



Pxy = (0.23,0.34,0.45,0.4), 
<jy = (0.3,0.25,0.35,0.5) 



Use d = (0.01,0.01,0.01, 
following solution: 



as the initial guess. The above algorithm then produces the 



D 



-0.06108 








0.2718 








-0.1145 










1 



R 



-0.1180 0.2300 

0.6466 -0.9388 

-0.5047 0.1955 

-0.5597 0.1657 



0.6490 -0.0990 

-0.7556 0.2306 

0.0117 -0.7905 

0.0880 0.5587 



Accordingly, in the Eq.CT p = (0.0678,0.0096,0.0062,0.0508) and b = -0.1446. The 
total time of calculation is 0.6 sec on a 32 bit PC with 3.0 Ghz single core CPU. 



5 Numerical method 



To solve the Eq. fll5p in general a iV-th dimensional variant of Strang's splitting 19j can be 
used which is 0(Ar 2 ). For linear operators this can be done by first formally solving the 
Eq. fTl5l) in the form 



$(r + Ar) = e AT ^ £l $(r) 



and then applying a generalized BCH formula 

At p At p A, 
g 2 *~0g 2 1 ...e 2 



For non-linear operators the situation is more delicate. However, as shown in [15] the pre- 
vious formal representation of the solution keeps to be valid in the non-linear well. 
Therefore, we can represent the previous equation as 



N 



8=1 
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and use the Strang's splitting. Explicitly this means that at each time step we have to solve 
a system of three equations 



Id 2 1 $*' 2 



1 



1 d 2 



1 $ 



2 u dy 2 



-fen ' yo 



G [0, Ar/2], 
G [0, Ar] 
G [0, Ar/2] 



(16) 



with the initial conditions for the first equation in Eq. ffl6|) : $*(0) = $(r), for the second 
one: $**(0) = $*(r + Ar/2), and for the last one: $(0) = $**(r + Ar). The final solution 
after this step is $(r + Ar) = $***(r + Ar/2). 

Since our terminal condition is of a rather complicated form given in the Eq. (|12|) . all equa- 
tions in Eq. (fl6l) can not be solved analytically, despite they do can be solved in quadratures. 
Indeed, the second equation is a iV-dimensional heat equation which admits an efficient nu- 
merical solution by using Fast Gauss Transform (FGT) since the Green's function is this 
case is a N- dimensional Gaussian. The remaining equations by change of variables known 
as Cole-Hopf transformation 



= a 6, $ = $i-o>o/»o) 



also reduces to the heat equation 



2 



■j/o.yo 



Therefore, they also can be solved by using FGT. 

Since we don't assume iV to be high, computation of the low-dimensional FGT doesn't 
face any difficulties if we use a powerful algorithm knows as Improved Fast Gauss Transform 
(IFGT) |21j |. Consider first a one-dimensional heat equation. Its solution can be represented 
as a convolution of the initial condition with the Green's function (which in this case is the 
Gaussian kernel). Suppose that the discretized space variable y is defined at M state nodes 
(source nodes). If we need to obtain the solution just at one fixed value of y , then we have 
one target point in space. However, according to the nature of the splitting algorithm we 
must solve similar problems at every splitting step (at given time we have 3 steps), and 
at every time step (the number of time steps J is determined as J = T/Ar). Therefore, 
to re-apply IFGT we need to use our target points as the initial points at the next step. 
Therefore, the number of the target points is also M. Then the total complexity of IFGT is 
0(2M). 

For d- dimensional problem the number of source and target points is M d . The complexity 
of IFGT is 0(2M d p(d) where f(d,p) is a polynomial function of d and the number of terms 
in (i-variate Taylor expansion truncated after order p — 1. To compare with finite-difference 
algorithms that usually are of the second order in space, consider an example with p=4 which 
provides a third order approximation. Thus, the total complexity of one step in time using 
Strang's splitting is 2M 5 (/(5,p) + 2/(1, p)). As shown in Q, e.g. /(5,4) = 56, /(1,4) = 
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4. Therefore, the complexity of the five- dimensional IFGT with M 5 source points and M 5 
target points is about 128M 5 . 

This could be compared with an analogous complexity of the finite difference method 
used to solve a d-dimensional heat equation at the space grid of M d nodes. Since all 
one-dimensional diffusion operators commute, this problem is reduced to five sequential 
one-dimensional problems. Every such a problem has the remaining M d ~ x states in other 
directions as dummy parameters, which means that this problem has to be solved M d ~ x 
times for every unique set of the dummy parameters. Also suppose we solve every problem 
with k steps in time (k = 6/ Ad)). Then the total complexity of the method is O(M) (the 
complexity of the FD one- dimensional solver for the heat equation, usually is about 6M) 
times k (the number of steps in time), times M d ~ x the number of the dummy variables) 
times d - the number of split tasks) which is 6kdM d . For d — 5 this gives 30/cM 5 . Therefore, 
at k > 4 IMGT is faster H At the same time the IMGT local error is essentially lower. That, 
as we mentioned, is because the standard schemes use the second order approximation in 
space 0, while the IFGT accuracy is defined by the number p, and is substantially higher. 

Accordingly, doing J steps in time results in the total complexity of the IFGT method 
to be 2JM N [f(N,p) + 2/(1, p)] = 0(JM N ). The proposed algorithm preserves the second 
order of approximation in time. 



6 Calibration 

To make this model practical one has to clearly understand how to calibrate the model to 
the market data. Two problems have to be discussed in this context. 

First we need to calibrate the risk-aversion parameter 7. Though this parameters may be 
specific to each investor, we may want to calibrate the risk aversion value to a "representa- 
tive" investor implied by the market. This problem was considered in |2J within a stochastic 
volatility model with a positive non-Gaussian Ornstein-Uhlenbeck process. Similar to our 
setup, the authors price options using the utility indifference with an exponential utility. 
The model is calibrated to historical returns, and the implied risk aversion is found by nu- 
merically inverting the indifference pricing equation given observed option prices. Certainly, 
in this case the risk aversion is a function of T and K, e.g. 7 = 7(T, K). 

An immediate problem with this approach is that when asset Yq is illiquid, it is hard to 
build the implied distribution of returns from the historical data, or to calibrate parameters 
of stochastic volatility for this asset. Therefore, in [2| liquid stocks (namely, MSFT and 
Volvo) were investigated. The initial intuition of the authors was that since the stochastic 
volatility model explains the observed market returns rather well, the implied risk aversion 
has to be almost flat with respect to T and K of the options. Contrary to this intuition, 
it was found that implied risk aversion exhibits a smile behavior for short dated options, 

5 Note that k=A is too small for any FD scheme to eliminate some additional errors produced by discon- 
tinuity in the first derivative of the payoff function. 

6 This produces a tri-diagonal matrix, and the total complexity of the solver is about 6M. Better approx- 
imations, e.g. using Pade schemes, lead to banded matrices, therefore the total complexity, while still linear, 
grows significantly (see, e.g. [l3j). 
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which was interpreted as issuers' fear of a market crash (in the case of the issuance of a put 
option). In particular, for Volvo, using call option bid/ask prices from December 30, 2005, it 
was found that risk aversion 7 varies from 0.1 to 0.01. It decreases when maturity increases 
from 1 month to 1 year, and also increases when K grows from 280 to 460. For puts the 
opposite is true, and the range of 7 is from 0.3 to 0. Similar behavior was observed for 
Microsoft, but in this case 7 reaches 10 for puts at K = 15 and T—l month. These results 
give an idea of a range of the implied risk aversion parameters. However, it doesn't address 
the above question of how to apply this approach to an illiquid asset. 

For some asset classes there sometimes exist other ways to imply the market value of 
7. For instance, for FX this problem is considered in 18j. An essential property of the FX 
market is the existence of cross- currency rules. For simple models of the underlyings (such as 
e.g. Geometric Brownian Motion, which is also our setting as well), this allows one to express 



7 in the explicit form via parameters of the domestic and foreign assets (see [18j|, Eq. 7.3.11). 
An example which uses monthly data for USD and GBP between December 31, 1985 and 
August 31, 2005, and DJI and FTSE as market representatives, gives an estimation 7=4.17. 
This implies that the choice 7=1 with a logarithmic utility function which is frequently used 
in the literature might not be very realistic. Note that a closed form expression for 7 is 
obtained in 18] for stochastic interest rates. 

Another challenge closely related to the first problem consists of the fact that for the 
illiquid asset Yq, it is hard to find its correlation with the potential candidates to be the 
proxy assets, Y 1 ,...,Y N , essentially almost by definition, as an illiquid asset typically does 
not move enough to measure its correlation with other assets. One way to proceed in such 
case is to use other, liquid assets from the same economic sector as Yq, as "correlation 
proxies" , as a way to roughly calibrate correlation parameters of Yq and our liquid proxies, 
which are the inputs in our framework. Note that in order to serve as a credible "correlation 
proxy" for Yq, another (liquid) proxy Y Q ' is expected to be similar to Yq, e.g. they should 
have similar credit ratings, credit default swap (CDS) spreads, expected default frequency 
(EDF) etc. 

The following differences of an illiquid asset from its liquid counterpart is discussed in 
First, an illiquid asset Yq can only be rebalanced at infrequent, stochastic intervals. When 
a trading opportunity arrives, the investor is able to rebalance her holdings of the illiquid 
asset. Furthermore, an illiquid asset is an asset that is not traded in a centralized exchange. 
In this case, investors who are willing to trade in this asset need to search for a counterparty. 
Such search process might be time-consuming, since in many cases the number of market 
participants with the required expertise, capital, and interest in these illiquid assets could 
be small. Examples of such illiquid assets are hedge funds, venture capital, private equity, 
structured credit, and real estate. Some of these assets are traded in OTC markets, but in 
others investors need to search directly for a counterparty in order to rebalance a position. 

The second way in which the illiquid asset differs from the liquid assets is that it cannot 
be pledged as collateral. Investors can issue non-state contingent debt by taking a short 
position in the riskless bond, but they cannot issue risky debt using the illiquid asset as 
collateral. If investors were allowed to do so, they could convert the illiquid asset into liquid 
wealth, and thus would implicitly circumvent the illiquidity friction. 
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This analysis means that the correlation between the illiquid asset Y and other proxy 
assets, at least in principal, can be computed from historical data, referring either Y (or its 
"correlation proxy" asset Yq). However, this is a delicate issue since the historical times series 
for Yq are recorded with time periods demonstrating kind of stochastic behavior. From this 
prospective an extended Kalman filter is a proper tool to work with the sparse, irregular time 
series. For more detail, see, e.g. {gj. Another prominent approach is a spectral estimation of 
a non-stationary time series sampled with missing data. The time series could be modeled 
as a locally stationary wavelet process, and its realization is assumed to feature missing 
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7 Conclusions 

In this paper we proposed a framework for pricing derivatives written on illiquid asset using a 
mixed dynamic-static hedging in a proxy index and N proxy options. While in this paper we 
apply our framework to an incomplete market version of the credit-equity Merton's model, 
the same approach can be used for other asset classes (equity, commodity, FX, etc.), e.g. for 
pricing and hedging options with illiquid strikes or illiquid exotic options, executive stock 
options etc. 

An efficient numerical algorithm is proposed which combines several changes of indepen- 
dent variables at the first step and Strang's splitting at the second step. 

A linear change of variables to new factorized (adiabatic) variables transforms the HJB 
equation for our model into a (N + l)-dimensional heat equation with an extra non-linear 
term. This term is proportional to $^ /$, i.e. it contains only the first derivative of the 
dependent variable $ wrt the first independent variable yo- This in contrast to the orig- 
inal HJB equation that has mixed derivatives, drifts and the non-linear term of the form 
(Ei=o ®yo) 2 /®- We propose an efficient numerical algorithm to compute coefficients of this 
linear transform. Some peculiarities of the algorithm are discussed. In particular, in the 
case of strong correlations between the illiquid asset Y and other proxy assets Yx, Y/v, the 
diagonal matrix D which we have to compute could be stiff. In this case, computation of 
eigenvectors of a non-symmetric matrix DA could require special methods (preconditioners) 
to preserve accuracy of computations. 

At the next step this new HJB equation in new variables is solved numerically using 
Strang's splitting. We show that this problem reduces to the solution of one iV-dimensional 
and two one- dimensional heat equations. Furthermore, we propose to use the Improved Fast 
Gauss Transform to decrease the total complexity of the method. We demonstrate that this 
complexity is 2JM N [f(N,p) + 2/(1, p)], where J is the number of steps in time, M is the 



number of grid points in S, and function f(m,ri) is defined in [21[. This algorithm is of the 
second order of approximation in time and of the p — 1 order of approximation in space. We 
also compare this with the finite-difference algorithm and show that our proposed algorithm 
produces less error and is more efficient in performance. 

In this paper for all assets we used a GBM model with time-dependent drifts and constant 
volatilities <7j. But this approach can also be generalized when volatilities <7j = <Ji{t) are 
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functions of time. This case will be discussed elsewhere. 
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A Proof of Proposition 14.1 

Proposition 14.11 claims that given a diagonal matrix D and a symmetric real matrix A, and 
matrix R of eigenvectors of DA, e.g. R~ X DAR = A, where A is a diagonal matrix with 
eigenvalues at the diagonal, it follows that R T AR is a diagonal matrix. 

Proof 2 By definition DA = RAR' 1 . Multiply both sides of this expression by D~ x l 2 from 
the left and by D 1 ! 2 from the right to obtain 

D^ 2 AD 1 ' 2 = D- l ' 2 RkR- l D 1 ' 2 

Introducing a diagonal matrix E such that EE -1 = I - a unit matrix, this can also be rewritten 
as 

D V I 2 AD 1 I 2 = D-^REE-'AR-'D 1 / 2 = D'^REAE^R^D 1 ' 2 . (17) 

Matrix D X I 2 AD X I 2 is a symmetric complex matrix, therefore it can be decomposed using its 
eigenvectors R and eigenvalues A which coincide with that of the matrix DA. 

D l/2 AD l/2 = RAR- 1 (18) 

Comparing the Eq. ( fj?| j and Eq. U~R) we see that eigenvectors R and R are connected by the 
map 

R = D^RE- 1 (19) 

Using this map and taking into account that D and E are diagonal matrices we can transform 
the matrix R T AR as follows 

R T AR = [D l l 2 RE- l ] T AD l l 2 RE- 1 = E~ l R T (D l / 2 ) T AD 1 ' 2 RE~ l 
= E- l R~ l D l ' 2 AD l l 2 RE- 1 = E^AE' 1 . 

Here we used the fact that the matrix D 1 ^ 2 is diagonal; D 1 ^ 2 AD 1 ! 2 is a complex symmetric 
matrix, therefore its eigenvectors R are orthogonal and R T = R' 1 . 

The last step of the proof is to recognize that since matrices A and E are diagonal, the 
product S _1 AS _1 is a diagonal matrix as well. 

Note, that matrix E is not an arbitrary matrix. It is determined by the Eq. / fT^j) and is 

E = R^D^R 

Accordingly, 

R T AR = R^D-^RAR^D-^R 
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